////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//---------------------------- Ebrahim Foulaadvand, 17 July 2013 -----------------                        //
//                                                                                                        //
//  The routine "Traffic" solves the 1D traffic equation using the FTCS, Lax and Lax-Wendroff methods.    //
//  Periodic boundary condition is used. The system length L is set to one. Initial wave shape is a step  //                                                                                               //
//  function.                                                                                             //         //
//                                                                                                        //                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>

using namespace std;

main()
{

double tau=0.001,h,L=1,x,vmax=1,rhomax=1;

int N=200,i,n,T=3000;


ofstream file0 ("initial profile n=0.plt");  //output file for the profile at timestep n=0.
ofstream file1 ("profile FTCS n=10.plt");    //output file for the profile at timestep n=100.
ofstream file2 ("profile FTCS n=20.plt");    //output file for the profile at timestep n=500.
ofstream file3 ("profile FTCS n=30.plt");    //output file for the profile at timestep n=500.
ofstream file4 ("profile Lax n=100.plt");     //output file for the profile at timestep n=500.
ofstream file5 ("profile Lax n=200.plt");    //output file for the profile at timestep n=500.
ofstream file6 ("profile Lax n=300.plt");    //output file for the profile at timestep n=500.
ofstream file7 ("profile LW n=1000.plt");     //output file for the profile at timestep n=500.
ofstream file8 ("profile LW n=2000.plt");    //output file for the profile at timestep n=500.
ofstream file9 ("profile LW n=3000.plt");   //output file for the profile at timestep n=500.
ofstream file10("profile analytic n=100.plt");
ofstream file11("profile analytic n=200.plt");
ofstream file12("profile analytic n=300.plt");



vector<double> rhoanal(N+1,0),rho(N+1,0),rhonew(N+1,0),rhoL(N+1,0),rhoLnew(N+1,0),rhoLW(N+1,0),rhoLWnew(N+1,0),
               rhohalfminus(N+1,0),rhohalfplus(N+1,0),chalfminus(N+1,0),chalfplus(N+1,0),J(N+1,0),
               JL(N+1,0),JLW(N+1,0);

 // Arrays "rho" and "rhonew" store the current and updated values of the density rho(x,t) at grid points.
 // "L" stands for Lax and "LW" stands for Lax-Wendroff.

h=L/N; cout<<"h= "<<h<<endl; cout<<"tau= "<<tau<<endl; getch();


for(i=0;i<=N;i++){ rho[i]=0; rhoL[i]=0; rhoLW[i]=0;

x=-0.5*L+i*h; if (i>=0.25*N && i<=0.5*N) { rho[i]=0.9*rhomax; rhoL[i]=0.9*rhomax; rhoLW[i]=0.9*rhomax; }

file0<<(-0.5*L+i*h)<<"  "<<rho[i]<<endl;




}

//storing the initial condition data in the data file.


for (n=1;n<=T;n++){ //cout<<"n= "<<n<<endl; //getch();


// -------------------       FTCS method -------------------------------------
/*
for(i=0;i<=N;i++){ J[i]=vmax*rho[i]*(1-(rho[i]/rhomax)); //cout<<"J["<<i<<"]="<<J[i]<<endl;getch();

}

rhonew[0]=rho[0]-(tau/(2*h))*(J[1]-J[N]); //cout<<"rhonew[0]="<<rhonew[0]<<endl;getch(); // periodic condition has been used: u[-1] is replaced with u[N]

for(i=1;i<N;i++){ rhonew[i]=rho[i]-(tau/(2*h))*(J[i+1]-J[i-1]); //cout<<"rhonew[i]="<<rhonew[i]<<endl;getch();

 }


rhonew[N]=rho[N]-(tau/(2*h))*(J[0]-J[N-1]); //cout<<"rhonew[N]="<<rhonew[N]<<endl;getch(); // periodic condition has been used: u[N+1] is replaced with u[0]


if (n==100)  for(i=0;i<=N;i++){ file1<<(-0.5*L+i*h)<<"  "<<rhonew[i]<<endl; } //storing in the data file.

if (n==200)  for(i=0;i<=N;i++){ file2<<(-0.5*L+i*h)<<"  "<<rhonew[i]<<endl; } //storing in the data file.

if (n==300)  for(i=0;i<=N;i++){ file3<<(-0.5*L+i*h)<<"  "<<rhonew[i]<<endl; } //storing in the data file.
*/

//-----------------------------     Lax Method      ------------------------------------

for(i=0;i<=N;i++){ JL[i]=vmax*rhoL[i]*(1-(rhoL[i]/rhomax)); }

rhoLnew[0]=0.5*(rhoL[1]+rhoL[N])-(tau/(2*h))*(JL[1]-JL[N]); // periodic condition has been used: u[-1] is replaced with u[N]

for(i=1;i<N;i++){ rhoLnew[i]=0.5*(rhoL[i+1]+rhoL[i-1])-(tau/(2*h))*(JL[i+1]-JL[i-1]); //cout<<"rhoLnew["<<i<<"]="<<rhoLnew[i]<<endl;getch();

}

rhoLnew[N]=0.5*(rhoL[0]+rhoL[N-1])-(tau/(2*h))*(JL[0]-JL[N-1]); // periodic condition has been used: u[N+1] is replaced with u[0]


if (n==100){

for(i=0;i<=N;i++){ x=-0.5*L+i*h; file4<<(-0.5*L+i*h)<<"  "<<rhoLnew[i]<<endl;  //storing in the data file.

if (x<=-vmax*tau*n) rhoanal[i]=0.9*rhomax;

if (x>-vmax*tau*n && x<vmax*tau*n) rhoanal[i]=0.5*0.9*rhomax*( 1-(x/(vmax*n*tau)) );

if (x>=vmax*tau*n) rhoanal[i]=0.0;

file10<<(-0.5*L+i*h)<<"  "<<rhoanal[i]<<endl;

}
} //storing in the data file.


if (n==200){

for(i=0;i<=N;i++){ x=-0.5*L+i*h; file5<<(-0.5*L+i*h)<<"  "<<rhoLnew[i]<<endl;

if (x<=-vmax*tau*n) rhoanal[i]=0.9*rhomax;

if (x>-vmax*tau*n && x<vmax*tau*n) rhoanal[i]=0.5*0.9*rhomax*( 1-(x/(vmax*n*tau)) );

if (x>=vmax*tau*n) rhoanal[i]=0.0;

file11<<(-0.5*L+i*h)<<"  "<<rhoanal[i]<<endl;

}
} //storing in the data file.


if (n==300){

for(i=0;i<=N;i++){ x=-0.5*L+i*h; file5<<(-0.5*L+i*h)<<"  "<<rhoLnew[i]<<endl;

if (x<=-vmax*tau*n) rhoanal[i]=0.9*rhomax;

if (x>-vmax*tau*n && x<vmax*tau*n) rhoanal[i]=0.5*0.9*rhomax*( 1-(x/(vmax*n*tau)) );

if (x>=vmax*tau*n) rhoanal[i]=0.0;

file12<<(-0.5*L+i*h)<<"  "<<rhoanal[i]<<endl;

}
} //storing in the data file.



//-----------------------------     Lax-Wendroff Method      ------------------------------------

for(i=0;i<=N;i++){ JLW[i]=vmax*rhoLW[i]*(1-(rhoLW[i]/rhomax)); }

rhohalfminus[0]=0.5*(rhoLW[N]+rhoLW[0]); rhohalfplus[0]=0.5*(rhoLW[1]+rhoLW[0]);

for(i=1;i<N;i++){ rhohalfminus[i]=0.5*(rhoLW[i-1]+rhoLW[i]); rhohalfplus[i]=0.5*(rhoLW[i+1]+rhoLW[i]); }

rhohalfminus[N]=0.5*(rhoLW[N-1]+rhoLW[N]); rhohalfplus[N]=0.5*(rhoLW[0]+rhoLW[N]);

//for(i=0;i<=N;i++){cout<<"rhohalfminus["<<i<<"]= "<<rhohalfminus[i]<<endl;  }

//for(i=0;i<=N;i++){cout<<"rhohalfplus["<<i<<"]= "<<rhohalfplus[i]<<endl;  }


for(i=0;i<=N;i++){ chalfminus[i]=vmax*(1-(2/rhomax)*rhohalfminus[i]); chalfplus[i]=vmax*(1-(2/rhomax)*rhohalfplus[i]); }

//for(i=0;i<=N;i++){cout<<"chalfminus["<<i<<"]= "<<chalfminus[i]<<endl;  }

//for(i=0;i<=N;i++){cout<<"chalfplus["<<i<<"]= "<<chalfplus[i]<<endl;  }




rhoLWnew[0]=rhoLW[0]-(tau/(2*h))*(JLW[1]-JLW[N]) + (tau*tau/(2*h*h))*( chalfplus[0]*(JLW[1]-JLW[0]) -
chalfminus[0]*(JLW[0]-JLW[N]) );


for(i=1;i<N;i++){

rhoLWnew[i]=rhoLW[i]-(tau/(2*h))*(JLW[i+1]-JLW[i-1]) + 0.5*(tau/h)*(tau/h)*( chalfplus[i]*(JLW[i+1]-JLW[i]) -
chalfminus[i]*(JLW[i]-JLW[i-1]) ); }


rhoLWnew[N]=rhoLW[N]-(tau/(2*h))*(JLW[0]-JLW[N-1]) + (tau*tau/(2*h*h))*( chalfplus[N]*(JLW[0]-JLW[N]) -
chalfminus[N]*(JLW[N]-JLW[N-1]) );


if (n==1000)  for(i=0;i<=N;i++){ file7<<(-0.5*L+i*h)<<"  "<<rhoLWnew[i]<<endl; } //storing in the data file.

if (n==2000)  for(i=0;i<=N;i++){ file8<<(-0.5*L+i*h)<<"  "<<rhoLWnew[i]<<endl; } //storing in the data file.

if (n==3000)  for(i=0;i<=N;i++){ file9<<(-0.5*L+i*h)<<"  "<<rhoLWnew[i]<<endl; } //storing in the data file.



for(i=0;i<=N;i++){ rho[i]=rhonew[i]; rhoL[i]=rhoLnew[i];  rhoLW[i]=rhoLWnew[i];  }


}


cout<<"end." <<endl; getch();

}